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ABSTRACT 

Extensive study on the numerical simulation of the vortical flow over 
a double-delta wing is carried out using the “thin-layer” Navier-Stokes and 
Euler equations. Two important flow characteristics, vortex interaction 
and vortex breakdown, are successfully simulated. Grid resolution is one 
of the most important factors associated with the vortex problem. Compu- 
tations were performed on a series of grids with various levels of refinement, 
coarse, medium and fine. Computations using either the coarse or medium 
grids fail to capture the proper physical phenomena. The computed result 
using a fine grid shows flow unsteadiness once the vortex breakdown takes 
place. The Cl - a characteristics are well predicted up to the breakdown 
angle of attack for all the grid distributions. The Euler solutions show 
fairly good agreement with the experiment on the Cl - a characteristics. 
However, other aspects of the solution at each angle of attack, such as the 
locus of the leading-edge separation vortex, are not consistent with the 
experiment. Even for the fine-grid Navier-Stokes computations, further 
grid resolution is required to obtain good quantitative agreement with the 
experiment. 

1. INTRODUCTION 

The vortical flow associated with leading-edge separation is believed 
to be essentially convection-dominated flow ( in other words, rotational in- 
viscid flow ), and thus, methods which describe flow separation and inviscid 
rotational flow are feasible for this purpose. Recently, many simulations of 
such flow fields using three-dimensional Euler and/or Navier-Stokes equa- 
tions have been reported (see Ref. 1). The computed results indicate that 
these methods are capable of predicting a variety of vortical flow fields. 

The Euler equations, which cannot describe flow separation math- 
ematically, give us reasonable results because the flow separates at the 
leading edge owing to the numerical dissipation. It is claimed that these 
solutions are essentially independent of the choice of dissipation, at least 
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for sharp leading-edge geometries [2]. One defect of the Euler equation ap- 
proach is that it cannot capture a secondary vortex near the leading edge. 
On the other hand, the Navier-Stokes equations have the capability to 
describe all the separation vortices. Most of computations have used the 
“thin-layer” Navier-Stokes equations; This is justified by the fact that the 
viscous effects are confined to a thin layer near the wall boundary and are 
dominated by the viscous terms associated with the strain rates normal to 
the wall, and the flow away from the body is essentially rotational inviscid. 
Viscous terms may not be negligible in the shear layer rolling up from the 
leading edge and in the core of the vortical flow. However, the contribu- 
tion of these terms to the basic structure of the flow field is believed to 
be small. Besides, the viscous terms are not properly evaluated in these 
regions, even by the Full Navier-Stokes equations, since the computational 
grid is usually not fine enough. Recent computations by Fujii et al. [3) 
and Thomas et al. [4] showed that vortex breakdown is simulated by the 
“thin-layer” Navier-Stokes equations. 

Although many computational results have been reported, questions 
remain to be answered about the reliability of the solutions. They are: 1) 
Is the solution quantitatively good? 2) Is the result independent of the 
number of the grid points to be used? 3) Is the result independent of the 
choice of numerical dissipation? 4) Does the Euler computation predict 
the vortex breakdown? 

In the present paper, three-dimensional Euler and Navier-Stokes 
computations for the nearly incompressible vortical flow fields over a 
double-delta wing are presented to help to answer these questions. The 
flow field over a double-delta wing is more complicated than that of a delta 
wing. There are separation vortices from both a strake leading edge and a 
wing leading edge. Over the wing region, these vortices interact with each 
other. The inboard vortex originally induced by the strake leading edge 
is relatively weak in the wing region, and thus it moves outboard by the 
induced velocity of the wing vortex and merging of the two vortices occurs. 
When the angle of attack is increased, vortex breakdown occurs near the 
wing trailing edge. To obtain numerical solutions for leading-edge vorti- 
cal flow fields, both a powerful computer and an efficient computational 
method are needed. Two supercomputers, the CRAY 2 at NASA Ames 
Research Center and the Amdahl 1200 at Amdahl Corporation, are used. 
The code based on the LU-ADI algorithm developed by one of the authors 
is appropriate for this purpose since extensive applications of this code for 
many other flow fields indicated reliability and efficiency. 

2. Governing Equations and 
Numerical Algorithm 

The partial differential equations in the generalized coordinate sys- 
tem, governing the three-dimensional flow of an unsteady, ideal gas can be 
written in conservation-law form: 
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d r Q + d^E + d„F + d<G = Re~ l d^S (1) 

The notations of each vector of Eq. (1) can be found in Ref. 3. In 
Eq. (1), the thin-layer approximation has been introduced. Euler equa- 
tions are obtained by neglecting the viscous terms in Eq. (1). The pressure, 
density, and velocity components are related to the energy by the following 
equation: 

p = (7 - l)(e - p(u 2 + v 2 + w 2 )/2) (2) 

In the following computations, flow is assumed to be laminar and no 
turbulence model is used. The metrics are evaluated using second-order 
central-difference formulas for interior points and three-point, one-sided 
formulas at the boundaries. Solution of Eqs. (1) is obtained by time 
integration. 

The numerical algorithm used here is the LU-ADI factorization 
method proposed by Fujii et al. [5] . Each ADI operator is decomposed 
into the product of lower and upper bidiagonal matrices by using a flux- 
vector-splitting technique and a diagonally dominant factorization. Cen- 
tral differencing is used in the right-hand side. In the solution process, 
an inversion in one direction consists of one forward sweep and one back- 
ward sweep of the scalar matrix. Thus, the LU-ADI algorithm requires 
little additional memory and is easily vectorized. The basic algorithm is 
first order in time and second order in space. For the convective terms 
in the right-hand side, fourth-order differencing is used except near the 
boundaries. Maintenance of the free stream is achieved by subtracting the 
free-stream fluxes from the governing equations. The artificial dissipation 
model developed in Ref. 5 is used in this study. 

3. RESULTS 

Most of the computations were carried out on the Amdahl 1200 
supercomputer with the maximum speed of 570 MFLOPS. Two of the 
fine-grid Navier-Stokes solutions showing vortex breakdown were obtained 
on the CRAY 2 supercomputer with the maximum speed of 1700 MFLOPS 
for four processors. One computation was carried out on both machines 
and the computed Cp distributions were within plotting accuracy. The 
code required 8.6 psec per grid point per iteration on the Amdahl 1200 
and 20.0 psec on a single processor of the CRAY 2. 

3.1 Grid Topology and Boundary Conditions 

The body geometry in the present study is a double-delta wing. 
The leading-edge swept angle is 80° at the strake and 60° at the wing. 
The thickness is 0.6% of the root chord and the leading edge is rounded. 
The geometry is the same as that of the experiment in Ref. 6, and is very 
similar to that of the experiment in Ref. 7 except for a slight difference of 
the wing thickness. A perspective view of the discretized region is shown in 
Fig. 1. The grid topology is H-type in the chordwise direction and O-type 
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in the spanwise direction. The computational grid is generated by the 
two-dimensional, hyperbolic, grid-generation method for each chordwise 
station. To accurately describe the body geometry, the surface grid is 
carefully distributed, especially in the apex, leading-edge, and trailing- 
edge regions. 

The grid extends approximately two root chords upstream, three 
chords downstream, and two and one-half chords above, below, and out- 
board of the wing. The location of the outer boundaries is not far from 
the body but is probably satisfactory. The influence of the outer bound- 
ary on the solutions needs further investigation. Free-stream values are 
specified along the upstream and circumferential boundaries and the pres- 
sure is fixed, and the extrapolation of other physical variables are used 
at the outflow boundary. Bilateral symmetry is imposed to reduce the 
computational domain. 

The grid consists of 119 points in the chordwise ( f ) direction, with 
14 points upstream of the body, 76 points over the body and 25 points in 
the wake; 101 points in the circumferential ( r) ) direction; and 71 points 
in the radial ( f ) direction. To see the effect of grid resolution, two other 
computational grids are used that consist of 61x63x41 points in the f, 77 , 
and f directions, respectively (medium grid), and 41x33x27 points (coarse 
grid). It should be noted that in both the medium and the coarse grids 
the computational surface grids lack representation of the geometry near 
the apex, the trailing edge, and certain areas of the leading edge because 
the total number of the grid points is limited. 

The computations reported in this paper utilize a Mach number of 
0.3 so that comparison with a low-speed experiment can be made. The 
Reynolds number is 1.3 x 10 6 based on the chord length at the root corre- 
sponding to the experiment[6]. 

3.2 Fine Grid Navier-Stokes Solutions 

The perspective view of the total pressure contour plots in several 
chordwise stations is shown in Fig. 2 at a=12°. Total pressure contour 
plots are frequently used in both experiments and computations because 
they are good indication of vortex patterns. At this angle of attack, two 
vortices, one from the strake and the other from the wing leading edge, are 
observed. The strength of the wing vortex increases with increasing dis- 
tance downstream of the strake-wing junction. As the strength of the wing 
vortex increases, the weaker strake vortex is drawn outward and downward 
owing to the velocity field induced by the wing vortex. At the downstream 
station the strake vortex is observed to merge with the wing vortex. The 
total pressure loss due to the secondary separation vortex is also seen in 
the two downstream stations shown in Fig. 2. The merging process is 
well simulated by the computation. However, the surface pressure plots 
(not shown here) indicate that the computed inner strake vortex is weaker 
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compared to the experiment^] . Downstream of the strake-wing junction, 
the strake vortex is no longer fed by vorticity shed from the leading edge, 
and so remains constant in strength, or may weaken because of viscous 
effects. The current grid resolution and the spatial accuracy of the com- 
putational scheme are not enough to capture this process accurately. In 
other words, the vortex tends to lose its strength rapidly because of the 
numerical dissipation. 


The perspective view of the total pressure contour plots at several 
chordwise stations is shown in Fig. 3 at a=30°. At this angle of attack, 
the merging of the strake and wing vortices occurs near the strake-wing 
junction. An abrupt change of the flow is observed at about 85 % of 
chord along the streamwise direction, where the size of the vortex core 
suddenly becomes large. Off-surface particle-path traces shown in Fig. 
4 demonstrate the breakdown of the primary vortex clearly. The flow 
pattern obtained in the experiment]?] for the same planform geometry is 
also presented. It should be noted that the experimental Reynolds number 
is lower and weak asymmetry of the flow is observed in the experiment. 
The vortex from the strake region is well-ordered and is very tightly coiled 
until about 80 or 85% of the chord, where an abrupt change of the flow 
is observed. In addition to the particles released from the strake and 
the wing leading edges, some particles are released from the trailing-edge 
region. Some of these particles move upstream to about 85% of chord and 
then move back downstream to the trailing-edge region. The existence of 
streamwise reverse flow is an evidence of the vortex breakdown. It should 
be noted that these reverse-flow particles continue to swirl along with the 
general motion of the vortex although the swirling is weak. It is also 
noticed that the strake vortex moves outward just behind the strake-wing 
junction because of the interaction with the wing vortex. 

The top view of the off-surface particle path traces at a=35° is 
shown in Fig. 5. The breakdown now takes place at about 60% chord 
location. The core flow from the strake leading edge has an abrupt kink 
there and forms a large spiral. A vertical spiral inside the breakdown region 
is observed at this location. This figure indicates that the core flow does 
not keep its original strong vorticity after the kink point. This may come 
from the resolution problem mentioned above. However, it is obvious that 
the a=35° case shows a different type of breakdown from the breakdown 
at a=30°. 


Careful examination of the flow field (see Ref. 3) revealed that the 
breakdown at a=30° can be classified to be a bubble-type breakdown. 
Even though there is no clear evidence, the result at o=35° seems to 
indicate the existence of the spiral- type breakdown. The structure of the 
vortex breakdown is not fully understood. In the experiment[8] on the 
delta wing, for instance, it is reported that both a spiral and a bubble 
type of breakdown were observed to transform one to the other at random. 
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The effect of turbulence is one factor always in issue. The computation 
assumes laminar flow, and the grid resolution is not enough to capture 
small eddies associated with the turbulent flow. However, global vortical 
flow structure in the absence of the breakdown must remain the same for 
both laminar and turbulent flows. Thus, the computed solutions may be 
useful to understand the flow structure, although the strict quantitative 
comparison with the experiment is very difficult. 

Figure 6 shows the Cl - a curve comparison for the present Navier- 
Stokes and Euler calculations and experimental data of Ref. 6. The 
Navier-Stokes result shows good agreement up to about 27° where vortex 
breakdown takes place. Careful examination of the Cl - a curve indicates 
nonlinear behavior of the flow field. The critical change of lift occurs at 
about 27° angle of attack with vortex breakdown. The exact angle of 
attack at which vortex breakdown occurs is sensitive to the initial condi- 
tion, probably because of a hysteresis effect. The Cl value in this figure 
is obtained with the free-stream initial condition. When started from the 
solution at a=30°, the computed value at a-28° becomes 0.912±0.005, 
which is less than the value in this figure. The Cl value at q= 27° becomes 
almost constant at the value 1.21 when started from the solution at q= 25°, 
and becomes 0.96610.004 when started from the solution at a=28°. The 
residual does not drop enough to obtain steady-state solutions for these 
cases as described below although the Cl values stay in a certain range 
even after 4000 to 8000 additional iterations. The outflow condition and 
the location of the outer boundary may explain this problem. More careful 
specification of the boundary conditions is necessary to discuss the quan- 
titative feature of the solutions. It should be remembered that a similar 
hysteresis of the Cl - ot curve is also observed in some experiments. The 
Euler solutions shown in Fig. 6 are discussed later. 

The convergence history for the <*=15° case is shown in Fig. 7 in 
terms of the Cl and the maximum residual. Since preliminary computation 
showed the flow unsteadiness at high angles of attack, the time step was 
set to be constant (At=0.01). The computation was continued up to 5000 
iterations for each case. At this angle of attack, the maximum residual 
drops more than three orders of magnitude within 2000 iterations and the 
Cl becomes almost constant. The Cl is not strictly constant because 
there is small unsteadiness in the flow field. However, the global flow 
field is settled and does not change. A similar tendency is obtained for 
every angle of attack up to the breakdown point. The computed C l shows 
clear unsteadiness above this angle of attack. However, the Cl stays in 
a certain range, as is shown in Fig. 6. It should be noted that all the 
results for an angle of attack above 30° are instantaneous rather than 
steady-state results that correspond to about 5000 iterations. At a=30° 
the unsteadiness is small and the global feature of the flow does not change. 
At a=35° flow inside the breakdown changes its feature in time. There 
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may be a concern about using a diagonal-form LU-AD1 solution method for 
unsteady computations since the algorithm is at best first-order-accurate 
and non-conservative in time. However, in the present study, quantitative 
study of the time-dependent response is not the primary interest. The 
main interest is associated with qualitative comparison of flow quantity 
amplitudes. In addition, the flow oscillations observed are slow and shock- 
free. The use of the diagonal form should be adequate. 

3.3 Fine-Grid Euler Solutions 

To evaluate the reliability of the Euler solution, several Euler com- 
putations were carried out under the same conditions. The code was the 
same as the Navier-Stokes computation, except that the viscous terms were 
deleted and the non-slip condition was not imposed. In addition, the same 
grid was used. The smoothing terms were also set to be similar. 

Figure 6 indicates that Euler results showed fairly good agreement 
with the experiment and may be in even better agreement than the Navier- 
Stokes solution throughout the breakdown region. The perspective view 
of the total pressure contour plots is shown in Fig. 8 at a=12°. The 
same contour level was used as in Fig. 2. The Euler solution shows only 
a very weak vortical flow (thus, Cl, is lower than in the Navier-Stokes 
solution), and the two vortices do not merge even at the trailing edge. 
The result at a=30° is shown in Fig. 9. Again, the contour level of Fig. 
3 was used. The Euler solution now shows a clear vortical flow but fails 
to predict vortex breakdown. However, at a=35° vortex breakdown takes 
place. Since solutions are sensitive to the initial condition when vortex 
breakdown takes place, it is difficult to draw conclusions. However, in the 
present case, up to the angle of attack where vortex breakdown takes place, 
the computed result using the Euler equations are not consistent with the 
experiment. 

3.4 The Grid and Smoothing Effect 

The Cl, - a curves for coarse-, medium- and fine-grid solutions are 
shown in Fig. 10. Although the results are almost the same up to the 
breakdown point, only the fine-grid solution indicates the change in slope 
due to the vortex merging. Although not shown here, two distinct vortices 
are not observed anywhere at 12° in the coarse- and the medium-grid 
solutions. The perspective views of the total pressure contour plots at 
30° are shown in Fig. 11 for the coarse-grid solution. The comparison of 
Fig. 11 and Fig. 3 indicates that the primary vortex becomes dissipative 
as the computational grid becomes coarse. In Fig. 3 a total pressure 
loss in the shear layer rolling up to the primary vortex is clearly seen. 
On the other hand, this is not observed in the coarse-grid solution. It is 
also noticed that the primary vortex is located outboard and downward 
when the computational grid is coarse. When the grid becomes finer, the 
resolution of the flow over the upper surface is improved, and the size of 
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the secondary separation vortex becomes larger. This secondary vortex 
moves the primary vortex inboard and upward. Thus, improved resolution 
of the flow changes the location of the primary vortex. The Cl - a curve 
shows only a small discrepancy because the primary vortex is located near 
the surface in the coarse-grid solution, even though the vortex is weak. 
The total pressure loss in the center of the primary vortex core is different 
for each grid result. In the conical Euler solutions, the total pressure loss 
and the location of the primary vortex are basically indifferent to the grid 
size and the amount of numerical smoothing[2]. This is not true in the 
Navier-Stokes solutions because the resolution of the flow changes the size 
of the secondary separation. The effect of resolution is more serious for the 
double-delta wing since the strake vortex in the wing region is relatively 
weak. Besides, vortex breakdown occurs only in the fine-grid solution. The 
particle-path traces indicate no reverse flows in the medium- and coarse- 
grid solutions. The computed Cl - a curve shows a good result, even 
for the coarse grid computation. Thus, coarse grid may be useful for the 
engineering purposes. However, it should be remembered that the solution 
may lack flow structure. 

To check the effect of the numerical smoothing terms, the coefficient 
of the smoothing terms was artificially increased to be three times larger. 
The result still displayed the existence of vortex breakdown and was close 
to the result in Fig. 3, except that the secondary separation vortex was 
smaller. This, along with the Euler solution in Fig. 9, indicates that the 
smoothing terms are not a key parameter in deciding flow-field structure. 
The viscous shear layers over the surface appear to be the most important 
factor in deciding the vortical flow structure. 

4. SUMMARY AND CONCLUSIONS 

Extensive study on the numerical simulation of the vortical flow over 
a double-delta wing was carried out using the “thin-layer” Navier-Stokes 
and Euler equations. Two important flow characteristics, vortex interac- 
tion and vortex breakdown, have been successfully simulated. The Cl - 
a characteristics were well predicted up to the breakdown angle of attack. 
The Euler solutions showed fairly good agreement with the experiment on 
the Cl - a characteristics. However, other aspects of the solution at each 
angle of attack such as the locus of the leading-edge separation vortex were 
not consistent with the experiment. Even for the Navier-Stokes computa- 
tions using a fine grid, further grid resolution is required to obtain good 
quantitative agreement with the experiment. Additional research into the 
careful specification of boundary conditions may also be required. 
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Fig. 2 Perspective view of the spanwise total pressure contour plots; 
Moo = 0.3, Re=1.3xl0 c \ a = 12.0°. 



h ig. 3 ] erspective view of the spanwise total pressure contour plots; 
Moo = 0.3, Re= 1.3x10°, n = 30.0°. 



Fig. 4 Computed off-surface particle-path pattern compared with the ex- 
perirnent[7]; M^ - 0.3, a = 30.0°, Re= ].3xl0 6 (Incompressible 
and Re=1.0xl0 5 for the experiment). 
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Fig. 6 Ci .vs. a characteristics; M ^ = 0.3, Re= 1.3x10°. 



ITERATION 

Fig. 7 Ci and maximum residual history for a = 15.0° ; A/oo - 0.3, 
Re=1.3xl0 6 . 


11 


Max Residual 





Pig. 8 Perspective view of the spanwise total pressure contour plots (E 
ler solution); - 0.3, a = 12.0°. 



F’ig- 9 Perspective view' of the spanwise total pressure contour plots (E 
ler solution); M ^ ----- 0.3, a = 30.0°. 
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